Gärtner et al. Plasmacytoid dendritic cells regulate tissue homeostasis of megakaryocytes

R Markdown

#Gärtner et al. Plasmacytoid dendritic cells regulate tissue homeostasis of megakaryocytes

library(SoupX)

## Warning: package 'SoupX' was built under R version 4.1.2

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(patchwork)
library(sctransform)

## Warning: package 'sctransform' was built under R version 4.1.2

library(ggplot2)
library(dittoSeq)
library(RColorBrewer)

## Warning: package 'RColorBrewer' was built under R version 4.1.2

library(data.table)

## 
## Attaching package: 'data.table'

## The following objects are masked from 'package:dplyr':
## 
##     between, first, last

library(magrittr)

## Warning: package 'magrittr' was built under R version 4.1.2

library(readxl)

## Warning: package 'readxl' was built under R version 4.1.2

library(reshape2)

## 
## Attaching package: 'reshape2'

## The following objects are masked from 'package:data.table':
## 
##     dcast, melt

library(enrichR)

## Warning: package 'enrichR' was built under R version 4.1.2

## Welcome to enrichR
## Checking connection ...

## Enrichr ... Connection is Live!
## FlyEnrichr ... Connection is Live!
## WormEnrichr ... Connection is Live!
## YeastEnrichr ... Connection is Live!
## FishEnrichr ... Connection is Live!
## OxEnrichr ... Connection is Live!

library(Seurat)

## Warning: package 'Seurat' was built under R version 4.1.2

## Attaching SeuratObject

library(clusterProfiler)

## Warning: package 'clusterProfiler' was built under R version 4.1.2

## 
## Registered S3 method overwritten by 'ggtree':
##   method      from 
##   identify.gg ggfun
## clusterProfiler v4.2.2  For help: https://yulab-smu.top/biomedical-knowledge-mining-book/
## 
## If you use clusterProfiler in published research, please cite:
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. The Innovation. 2021, 2(3):100141
## 
## Attaching package: 'clusterProfiler'
## 
## The following object is masked from 'package:stats':
## 
##     filter

library(stringr)

## Warning: package 'stringr' was built under R version 4.1.2

library(EnhancedVolcano)

## Loading required package: ggrepel

library(dittoSeq)
library(ComplexHeatmap)

## Loading required package: grid
## ========================================
## ComplexHeatmap version 2.10.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
## 
## If you use it in published research, please cite:
## Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
##   genomic data. Bioinformatics 2016.
## 
## The new InteractiveComplexHeatmap package can directly export static 
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
## 
## This message can be suppressed by:
##   suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================

library(RColorBrewer)
library(circlize)

## Warning: package 'circlize' was built under R version 4.1.2

## ========================================
## circlize version 0.4.15
## CRAN page: https://cran.r-project.org/package=circlize
## Github page: https://github.com/jokergoo/circlize
## Documentation: https://jokergoo.github.io/circlize_book/book/
## 
## If you use it in published research, please cite:
## Gu, Z. circlize implements and enhances circular visualization
##   in R. Bioinformatics 2014.
## 
## This message can be suppressed by:
##   suppressPackageStartupMessages(library(circlize))
## ========================================

MKP3HTO = readRDS("~/Downloads/final_MKP3HTO.RDS")
DefaultAssay(MKP3HTO) <- "RNA"

# Projecting singlet identities on UMAP  visualization
print(DimPlot(MKP3HTO, group.by = "HTO_classification"))

print(DimPlot(MKP3HTO, shuffle = T, seed = 1, group.by= "HTO_classification", split.by= "HTO_classification", ncol=3, raster=FALSE))

print(DimPlot(MKP3HTO, group.by= "seurat_clusters", raster=FALSE, label = TRUE))

MKP3HTO_markers <- FindAllMarkers(MKP3HTO, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25, test.use = "wilcox")

## Calculating cluster Mixed_progenitors

## Calculating cluster Late_MKP

## Calculating cluster MK-MEP

## Calculating cluster MK-MEP_cycling

## Calculating cluster GMP

## Calculating cluster Early_MKP

## Calculating cluster Immature_neutrophils

## Calculating cluster Basophils

## Calculating cluster NK

## Calculating cluster Erythroid

MKP3HTO_markers %>%
    group_by(cluster) %>%
    slice_max(n = 2, order_by = avg_log2FC)

## # A tibble: 20 × 7
## # Groups:   cluster [10]
##        p_val avg_log2FC pct.1 pct.2 p_val_adj cluster              gene     
##        <dbl>      <dbl> <dbl> <dbl>     <dbl> <fct>                <chr>    
##  1 0               1.39 0.894 0.494 0         Mixed_progenitors    AY036118 
##  2 0               1.13 0.996 0.976 0         Mixed_progenitors    Camk1d   
##  3 0               4.31 0.965 0.358 0         Late_MKP             Tsc22d1  
##  4 0               4.14 0.999 0.688 0         Late_MKP             Ctla2a   
##  5 4.45e-252       1.28 0.835 0.223 8.69e-248 MK-MEP               Apoe     
##  6 7.00e-178       1.18 0.987 0.715 1.37e-173 MK-MEP               Pdcd4    
##  7 8.48e-199       1.27 0.944 0.389 1.66e-194 MK-MEP_cycling       H2afx    
##  8 6.68e-180       1.25 0.898 0.383 1.30e-175 MK-MEP_cycling       Hist1h2ap
##  9 0               3.32 0.924 0.219 0         GMP                  Lgals1   
## 10 1.73e-246       3.13 0.9   0.414 3.38e-242 GMP                  Prtn3    
## 11 1.19e-176       2.58 0.915 0.326 2.33e-172 Early_MKP            Prkg1    
## 12 2.63e-157       2.42 0.994 0.551 5.13e-153 Early_MKP            Pbx1     
## 13 1.29e-137       7.16 0.918 0.373 2.51e-133 Immature_neutrophils S100a9   
## 14 7.94e-120       6.90 0.828 0.319 1.55e-115 Immature_neutrophils S100a8   
## 15 0               7.06 0.855 0.102 0         Basophils            Prss34   
## 16 0               6.15 0.897 0.09  0         Basophils            Mcpt8    
## 17 1.17e-255       4.90 0.453 0.024 2.29e-251 NK                   Ccl5     
## 18 0               3.43 0.608 0.009 0         NK                   Trbc2    
## 19 1.05e- 51       6.70 0.679 0.333 2.05e- 47 Erythroid            Hbb-bs   
## 20 1.43e-101       6.38 0.69  0.173 2.80e- 97 Erythroid            Hbb-bt

# plotting the marker genes
MKP3HTO_markers %>%
    group_by(cluster) %>%
    top_n(n = 10, wt = avg_log2FC) -> top10
print("Dotplot of cluster defining genes")

## [1] "Dotplot of cluster defining genes"

print(DotPlot(MKP3HTO, features = unique(top10$gene))+ RotatedAxis())

gene set enrichment analysis using the enrichR package for the top 50 marker genes of each cluster to better understand the composition of each cluster.

Idents(MKP3HTO) = MKP3HTO$seurat_clusters
new.cluster.ids <- c("Mixed_progenitors","Late_MKP","MK-MEP","MK-MEP_cycling","GMP","Early_MKP","Immature_neutrophils","Basophils","NK","Erythroid")
names(new.cluster.ids) <- levels(MKP3HTO)
MKP3HTO <- RenameIdents(MKP3HTO, new.cluster.ids)
MKP3HTO$Idents = Idents(MKP3HTO)

MKP3HTO_markers <- FindAllMarkers(MKP3HTO, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25, test.use = "wilcox")

## Calculating cluster Mixed_progenitors

## Calculating cluster Late_MKP

## Calculating cluster MK-MEP

## Calculating cluster MK-MEP_cycling

## Calculating cluster GMP

## Calculating cluster Early_MKP

## Calculating cluster Immature_neutrophils

## Calculating cluster Basophils

## Calculating cluster NK

## Calculating cluster Erythroid

MKP3HTO_markers %>%
    group_by(cluster) %>%
    slice_max(n = 2, order_by = avg_log2FC)

## # A tibble: 20 × 7
## # Groups:   cluster [10]
##        p_val avg_log2FC pct.1 pct.2 p_val_adj cluster              gene     
##        <dbl>      <dbl> <dbl> <dbl>     <dbl> <fct>                <chr>    
##  1 0               1.39 0.894 0.494 0         Mixed_progenitors    AY036118 
##  2 0               1.13 0.996 0.976 0         Mixed_progenitors    Camk1d   
##  3 0               4.31 0.965 0.358 0         Late_MKP             Tsc22d1  
##  4 0               4.14 0.999 0.688 0         Late_MKP             Ctla2a   
##  5 4.45e-252       1.28 0.835 0.223 8.69e-248 MK-MEP               Apoe     
##  6 7.00e-178       1.18 0.987 0.715 1.37e-173 MK-MEP               Pdcd4    
##  7 8.48e-199       1.27 0.944 0.389 1.66e-194 MK-MEP_cycling       H2afx    
##  8 6.68e-180       1.25 0.898 0.383 1.30e-175 MK-MEP_cycling       Hist1h2ap
##  9 0               3.32 0.924 0.219 0         GMP                  Lgals1   
## 10 1.73e-246       3.13 0.9   0.414 3.38e-242 GMP                  Prtn3    
## 11 1.19e-176       2.58 0.915 0.326 2.33e-172 Early_MKP            Prkg1    
## 12 2.63e-157       2.42 0.994 0.551 5.13e-153 Early_MKP            Pbx1     
## 13 1.29e-137       7.16 0.918 0.373 2.51e-133 Immature_neutrophils S100a9   
## 14 7.94e-120       6.90 0.828 0.319 1.55e-115 Immature_neutrophils S100a8   
## 15 0               7.06 0.855 0.102 0         Basophils            Prss34   
## 16 0               6.15 0.897 0.09  0         Basophils            Mcpt8    
## 17 1.17e-255       4.90 0.453 0.024 2.29e-251 NK                   Ccl5     
## 18 0               3.43 0.608 0.009 0         NK                   Trbc2    
## 19 1.05e- 51       6.70 0.679 0.333 2.05e- 47 Erythroid            Hbb-bs   
## 20 1.43e-101       6.38 0.69  0.173 2.80e- 97 Erythroid            Hbb-bt

MKP3HTO_markers %>%
    group_by(cluster) %>%
    top_n(n = 50, wt = avg_log2FC) -> top50

top50 = as.data.frame(top50)
conditions <- c("Mixed_progenitors","Late_MKP","MK-MEP","MK-MEP_cycling","GMP","Early_MKP","Immature_neutrophils","Basophils","NK","Erythroid")
for (ident in conditions){
enrich_genes = top50[top50$cluster == ident, ]$gene
    enrichment_results <- enrichR::enrichr(
      genes = enrich_genes,
      databases = "GO_Biological_Process_2023"
    )
    enrichment_results = as.data.frame(enrichment_results)
    enrichment_data <- data.frame(
  Term = enrichment_results$GO_Biological_Process_2023.Term,
  Overlap = enrichment_results$GO_Biological_Process_2023.Overlap,
  P.value = enrichment_results$GO_Biological_Process_2023.P.value,
  Adjusted.P.value = enrichment_results$GO_Biological_Process_2023.Adjusted.P.value,
  Odds.Ratio = enrichment_results$GO_Biological_Process_2023.Odds.Ratio,
  Combined.Score = enrichment_results$GO_Biological_Process_2023.Combined.Score,
  Genes = enrichment_results$GO_Biological_Process_2023.Genes
)
       enrich = plotEnrich(enrichment_data)
       print(ident)
       print(enrich)
}

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.

## Warning in plotEnrich(enrichment_data): There are duplicated trimmed names in
## the plot, consider increasing the 'numChar' setting.

## [1] "Mixed_progenitors"

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.
## [1] "Late_MKP"

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.
## [1] "MK-MEP"

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.
## [1] "MK-MEP_cycling"

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.

## Warning in plotEnrich(enrichment_data): There are duplicated trimmed names in
## the plot, consider increasing the 'numChar' setting.

## [1] "GMP"

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.
## [1] "Early_MKP"

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.
## [1] "Immature_neutrophils"

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.

## Warning in plotEnrich(enrichment_data): There are duplicated trimmed names in
## the plot, consider increasing the 'numChar' setting.

## [1] "Basophils"

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.

## Warning in plotEnrich(enrichment_data): There are duplicated trimmed names in
## the plot, consider increasing the 'numChar' setting.

## [1] "NK"

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.
## [1] "Erythroid"

getting proportion of cells per cluster belonging to each of the 4 conditions using the dittoBarPlot package.

print("getting proportion of cells per cluster for each condition")

## [1] "getting proportion of cells per cluster for each condition"

#BiocManager::install("dittoSeq")


dittoBarPlot(
    object = MKP3HTO,
    var = "CSclassification",
    group.by = "RNA_snn_res.0.2")

reading in the top 40 response to type I interferon genes from the bulk-seq reactome pathway analysis

DefaultAssay(MKP3HTO) <- "RNA"
## response to type I interferon genes as bulk_seq genes
bulk_seq = read.csv("~/Downloads/GOBP_RESPONSE_TO_TYPE_I_INTERFERON.csv", header = TRUE)
print("response to interferon genes")

## [1] "response to interferon genes"

head(bulk_seq)

##    NAME  PROBE
## 1 row_0   TYK2
## 2 row_1  HLA-A
## 3 row_2  PTPN6
## 4 row_3 CACTIN
## 5 row_4   BST2
## 6 row_5 YTHDF2
##                                                         DESCRIPTION.br..from.dataset.
## 1                               tyrosine kinase 2 [Source:HGNC Symbol;Acc:HGNC:12440]
## 2     major histocompatibility complex, class I, A [Source:HGNC Symbol;Acc:HGNC:4931]
## 3 protein tyrosine phosphatase non-receptor type 6 [Source:HGNC Symbol;Acc:HGNC:9658]
## 4           cactin, spliceosome C complex subunit [Source:HGNC Symbol;Acc:HGNC:29938]
## 5               bone marrow stromal cell antigen 2 [Source:HGNC Symbol;Acc:HGNC:1119]
## 6    YTH N6-methyladenosine RNA binding protein 2 [Source:HGNC Symbol;Acc:HGNC:31675]
##   GENE.SYMBOL GENE_TITLE RANK.IN.GENE.LIST RANK.METRIC.SCORE RUNNING.ES
## 1        null       null                54          2.885287 0.04354477
## 2        null       null               133          2.323923 0.07665094
## 3        null       null               194          2.121879 0.10751816
## 4        null       null               234          2.037717 0.13822223
## 5        null       null               255          1.982451 0.16911610
## 6        null       null               289          1.911551 0.19812344
##   CORE.ENRICHMENT
## 1             Yes
## 2             Yes
## 3             Yes
## 4             Yes
## 5             Yes
## 6             Yes

## getting the top 40 genes from response to type I interferon genes
bulk_seq = as.data.frame(bulk_seq)
rownames(bulk_seq) = bulk_seq$PROBE
bulk_seq_genes = head(rownames(bulk_seq), 40)
## converting the genes to mouse gene form - "Tyk2" "Hla-a"
bulk_seq_genes = tolower(bulk_seq_genes)
firstup = function(x) {
  substr(x, 1, 1) = toupper(substr(x, 1, 1))
  
  x
}
bulk_seq_genes = firstup(bulk_seq_genes)
print("response to interferon genes after processing")

## [1] "response to interferon genes after processing"

bulk_seq_genes

##  [1] "Tyk2"     "Hla-a"    "Ptpn6"    "Cactin"   "Bst2"     "Ythdf2"  
##  [7] "Hsp90ab1" "Irf9"     "Ythdf3"   "Ube2k"    "Rnasel"   "Setd2"   
## [13] "Shmt2"    "Cnot7"    "Irf1"     "Sp100"    "Irak1"    "Abce1"   
## [19] "Mul1"     "Egr1"     "Adar"     "Ifitm1"   "Rsad2"    "Jak1"    
## [25] "Irf5"     "Ifitm2"   "Ptpn2"    "Tbk1"     "Dcst1"    "Oasl"    
## [31] "Mavs"     "Xaf1"     "Trim6"    "Irf2"     "Ifna1"    "Ifnar2"  
## [37] "Ifitm3"   "Usp18"    "Cdc37"    "Ifna16"

### renaming the idents of seurat obj to conditions instead of clusters to calculate fold change between conditions for the genes
Idents(MKP3HTO) = MKP3HTO$CSclassification

projecting the response to type I interferon genes onto the single cell data using AddModuleScore function

print("projecting the response to type I interferon genes onto the single cell data")

## [1] "projecting the response to type I interferon genes onto the single cell data"

MKP3HTO <- AddModuleScore(MKP3HTO,
                    features = list(bulk_seq_genes), name = "bulk_seq_genes")

## Warning: The following features are not present in the object: Hla-a, Oasl,
## Ifna1, Ifna16, not searching for symbol synonyms

print(VlnPlot(MKP3HTO, features = "bulk_seq_genes1", group.by ="CSclassification"))

print(VlnPlot(MKP3HTO, features = "bulk_seq_genes1", split.by ="CSclassification"))

## The default behaviour of split.by has changed.
## Separate violin plots are now plotted side-by-side.
## To restore the old behaviour of a single split violin,
## set split.plot = TRUE.
##       
## This message will be shown once per session.

print(VlnPlot(MKP3HTO, features = "bulk_seq_genes1", group.by ="CSclassification", split.by = "RNA_snn_res.0.2"))

print(VlnPlot(MKP3HTO, features = "bulk_seq_genes1", split.by ="CSclassification", group.by = "RNA_snn_res.0.2"))

print(VlnPlot(MKP3HTO, features = "bulk_seq_genes1", group.by = "RNA_snn_res.0.2"))

identifying the DE genes from bulk-seq data by determining the cutoff as log FC > 0.25 and FDR < 0.05 per comparison. We then projected these bulk seq DE genes onto the single cell data for each comparison using the Module Score Function.

print("calculating DE genes from bulk_seq data as log FC > 0.25 and FDR < 0.05")

## [1] "calculating DE genes from bulk_seq data as log FC > 0.25 and FDR < 0.05"

bulk_de_genes = read.csv("~/Downloads/MKP_samples_log2FC_forVisha.csv", header = TRUE)
rownames(bulk_de_genes) = bulk_de_genes$Gene_name
# Create a subset of the data frame
genes_pdpdc_vs_cntrl  <- subset(bulk_de_genes, abs(log2FC_PD_DT_vs_Ctrl) > 0.25 & (FDR_PD_DT_vs_control) < 0.05)
genes_pdpdc_vs_cntrl = rownames(genes_pdpdc_vs_cntrl)

genes_pd_vs_pdpdc  <- subset(bulk_de_genes, abs(log2FC_PD_vs_PD_DT) > 0.25 & (FDR_PD_vs_PD_DT) < 0.05)
genes_pd_vs_pdpdc = rownames(genes_pd_vs_pdpdc)

genes_pd_vs_cntrl  <- subset(bulk_de_genes, abs(log2FC_PD_vs_Ctrl) > 0.25 & (FDR_PD_vs_control) < 0.05)
genes_pd_vs_cntrl = rownames(genes_pd_vs_cntrl)

print("projecting the bulk rna DE genes (platelet depletion vs Bl6) onto the single cell data")

## [1] "projecting the bulk rna DE genes (platelet depletion vs Bl6) onto the single cell data"

genes_pd_vs_cntrl_pos  <- subset(bulk_de_genes, (log2FC_PD_vs_Ctrl) > 0.25 & (FDR_PD_vs_control) < 0.05)

MKP3HTO <- AddModuleScore(MKP3HTO,
                    features = list(genes_pd_vs_cntrl_pos), name = "genes_pd_vs_cntrl_pos")

## Warning: The following features are not present in the object: c("Emp3", "Sf3b1", "Ccnk", "Idh3a", "Prtn3", "Ndufs8", "Ssrp1", "Hspa4", "Akap8", "Ctsg", "Ccnd2", "Smim14", "Clk4", "Capza1", "Fbxl14", "Nudcd2", "1110004F10Rik", "Wdr81", "Cnn2", "Hnrnpdl", "Mafg", "Naa38", "Dusp2", "Ivns1abp", "Dnmt1", "Mcm3", "Mrpl43", "Mbtd1", "1810021M19Rik", "Sdhb", "Mcpt8", "Mir17hg", "Hars", "Mrpl34", "Rfc2", "Smc3", "Atad2", "Prg2", "Cited2", "1110017D15Rik", "Ywhab", "Itgb3", "Hnrnpa1", "Psmb7", "Ccnl2", "Mrpl16", "Crebzf", "D2Wsu81e", "Copa", "Ctdspl2", "Pcbp2", "Nefh", 
## "Cops3", "Lcp1", "Sypl", "Dnajb11", "Noc4l", "Gm40508", "Anp32e", "Shfm1", "Nhp2", "LOC108167358", "Mpo", "Tspan31", "Echs1", "Aoc3", "Cct7", "Myh9", "Ccdc115", "Arpc5", "Glipr1", "C330027C09Rik", "Cyp2u1", "Canx", "Clcnka", "Pdap1", "Stxbp1", "Dctd", "Pcyox1l", "Gm15264", "Kdm6b", "Ccr2", "Ifitm1", "Sf3b2", "Fdxr", "Slc44a4", "Mxd1", "Eif1a", "Plaur", "Larp4b", "Eif4a3", "Sgol1", "Mir703", "Nfic", "Ddx6", "Cggbp1", "Chd4", "Bzw1", "Etohd2", "Trub2", "Tyk2", "Igbp1", "Gng11", "Tmem88", "Ube2r2", 
## "Idh1", "Btg1", "Nup88", "Nsun2", "Ctsl", "Gm40227", "Rgs1", "AI467606", "Svip", "Cct6a", "Snhg5", "Ubc", "Erp29", "Etf1", "Tsfm", "Dnajc15", "Psmd8", "0610012G03Rik", "Mtch1", "Sel1l", "Nip7", "2610001J05Rik", "Pdhb", "Rtcb", "Adsl", "Slc4a11", "Exosc8", "Tctn2", "Lmo4", "Atp2a2", "Psmb1", "Rad9b", "Vim", "Pdia3", "Trim28", "Sppl2a", "Clec12a", "Idh2", "Vimp", "Rab8b", "Ccdc163", "Fam71e2", "Ndufs6", "Smim7", "Usp1", "Cdk8", "Gm5812", "2610008E11Rik", "Trmt2a", "Mrfap1", "Gm4617", "Gm12992", "Snhg1", 
## "Nudt21", "Tdpx-ps1", "Cebpb", "Med13", "Mrpl20", "Mrpl53", "Laptm5", "Stk4", "Nsa2", "Pdia6", "Rrs1", "Ppp1r12a", "Apex1", "Cops8", "2700094K13Rik", "Tex261", "Cd47", "Tmem53", "Gnb2", "Trak1", "Itgal", "Tmem9", "Gm36368", "Helt", "Bcas2", "Spc25", "Tsnax", "Rabggtb", "Rsrc2", "Usp33", "Tufm", "2310009A05Rik", "LOC105247125", "Pam16", "Gm4737", "Fcgr3", "Cramp1l", "Gm42312", "Uap1", "5730422E09Rik", "B230206L02Rik", "Cdk1", "Gls", "Wdr83os", "Emc10", "Ppp6c", "Ccar1", "Gm36264", "Lyar", "Fam129a", 
## "Man2a1", "Timm17b", "1810026B05Rik", "Gm36876", "Xbp1", "Yy2", "Mien1", "Ube2b", "Fam107b", "Atp2b2", "Usf1", "Aurkaip1", "Prpf8", "Nup62", "Hmgn5", "Ndufb2", "Tmem14c", "Kpna2", "Mpeg1", "Lmna", "Trp53", "Bcl6", "Colgalt1", "Ndufb6", "Sap18b", "Ms4a3", "Rad51d", "Elane", "Slc17a9", "Cacybp", "Slc35b1", "Ubxn4", "Ddx3y", "H2afy", "Fus", "Nedd4", "Rexo2", "Stub1", "Metap2", "Gm31619", "Ddx26b", "Ube3c", "Cycs", "Atp6v1g1", "Cdv3", "Ube2k", "Tubb4b", "Gm41847", "4933404K13Rik", "Vbp1", "Aasdhppt", 
## "Cabp1", "Magt1", "Eif4g2", "Polr2k", "Bsn", "Rbmxl1", "Arhgef1", "Gm32673", "P2rx1", "Rnf5", "Ubald2", "Kif22", "Snrpd3", "Tjp2", "Ctsd", "Fbxl16", "LOC108168971", "Mlec", "Brd2", "Lage3", "Hmga1", "Cdk4", "Fam98b", "Pign", "Ier2", "Mrpl18", "Fdx1l", "Slc16a1", "Rplp0", "Psmb2", "Nktr", "Magoh", "Gm39923", "Eif2s2", "Ldb1", "Eif5b", "Pdcd5", "Ube2s", "Lsm4", "Osm", "Eif4g1", "Rgs2", "Ptprcap", "Malat1", "Pln", "Adprh", "Gstm1", "Rbm25", "Reep5", "Eef2", "Mc1r", "Tubb5", "Serpinb1a", "Wdr74", "Gtpbp2", 
## "Selenbp1", "Got2", "Timm13", "Itm2c", "Gm13394", "Irf2bp2", "Prcc", "Tnfaip2", "Ncl", "Tusc5", "Cdca3", "Acad11", "Slc25a3", "Ppp1cb", "Faap100", "Aars", "Junb", "Tnks2", "Id2", "Rbx1", "Srgn", "Ankle1", "Pea15a", "Gm15452", "Pltp", "Gm40017", "Gm3740", "Nedd8", "Ska1", "Psmg4", "5031434O11Rik", "LOC108169013", "Psmg2", "Srp19", "S100a9", "Gm36673", "Trem3", "Rps15a-ps4", "Srsf6", "Sf3a2", "Sh3glb2", "Csnk1g2", "E2f1", "Rasal3", "Rnaseh2c", "Tuba8", "Dynlt1b", "Nfix", "Mmadhc", "Rasgrp2", "Myb", 
## "Cenpl", "Brcc3", "Coasy", "Eea1", "Mrpl45", "Rnf180", "Gm14176", "Chchd7", "Stt3a", "Appbp2", "Gm20075", "E2f4", "Impdh2", "Stx3", "Mogs", "Skap2", "Prpsap2", "Gm10532", "Cnot1", "Eif3j1", "Gm6361", "1700123O20Rik", "Atp1b3", "Ndufs2", "Pdp2", "Psma6", "LOC108167318", "Jak3", "1110018N20Rik", "Mrpl57", "Tmem179b", "Pkn2", "Aamp", "Tcerg1", "Edem2", "Naa50", "1810022K09Rik", "Plpp5", "Rad21", "Rpl37a", "Rps2", "Zc3h7a", "Arih1", "Pcnp", "Tomm40", "Phf5a", "Tmed3", "Trmt1l", "Tm9sf1", "Atp5c1", "Eef1b2", 
## "Cldn13", "Plod3", "Cd164", "Bet1l", "Wdr59", "Csgalnact2", "Atr", "Emc6", "Ptp4a2", "Lsm7", "Gm6532", "Cops6", "Thoc2", "Mcm7", "Tyw5", "Ap2m1", "Gm10705", "Gtpbp4", "Eya2", "Tjp3", "Tdp2", "Ppib", "Acin1", "Mcrs1", "Gle1", "Rsl24d1", "F13a1", "Cab39", "Anp32b", "Rpl8", "Tbrg4", "Rps6", "Hmbs", "Nup85", "Nxf1", "Pdcd6", "Dbi", "Hmgcr", "Khsrp", "Abce1", "Cystm1", "Gm39351", "Tmem70", "Usp9x", "Fosl2", "Ier3", "Ireb2", "LOC108168421", "Neb", "Nol12", "Ubac1", "Ube2u", "Uevld", "Hn1", "Setd1b", "Irs2", 
## "Cox7a1", "Gm32534", "Actb", "Rpl7a", "Gm10444", "Lta4h", "Bola1", "Cul2", "Mex3c", "Pcsk4", "Hnrnpf", "Ap2s1", "Cbs", "Azin1", "Nmd3", "Fam120a", "Ube2a", "Inip", "Psma7", "Ncf2", "Eno1", "Rrnad1", "Rps11-ps2", "Uqcrh", "Atp5j", "Ccne2", "Eef1a1", "Mynn", "Tet1", "Mysm1", "Wsb1", "Syncrip", "Lamtor5", "Alkbh1", "Pigyl", "Gm32399", "4833439L19Rik", "Ahcy", "Arid4b", "Gm31194", "Irx5", "1110001J03Rik", "Hsp90ab1", "Sf3a1", "Ppm1g", "Qrich1", "Akna", "Klf6", "Pold1", "Tmem258", "Cinp", "Degs1", "Cwc15", 
## "Flna", "Gm39877", "Gm30718", "Glis1", "Polr2a", "Ptpn6", "Ranbp2", "Uba2", "Pabpc4", "A230103J11Rik", "Hsbp1", "Trim65", "Gm11914", "Osbp", "Csnk1g1", "Mrpl51", "Rbl1", "Zfp871", "Gdap10", "Kifc5b", "Zbtb7a", "Cep95", "Cuta", "Tctex1d2", "Rps12-ps10", "Ube2q1", "Pcm1", "Psmd4", "Capza2", "Gfi1", "Zfp825", "Gm30181", "Ces2g", "Gm5620", "Supt4a", "Snw1", "Emilin1", "Gnl3", "Mrps22", "Rpl13a", "Sh2b1", "Hspd1", "Gm41320", "Stard7", "Gm3386", "Kdelr2", "Olfr54", "Zbtb21", "Tsen54", "Acta2", "Sdha", 
## "Tmem170b", "5031434C07Rik", "Plagl2", "Nop16", "Atp5k", "Cks1b", "Cntrl", "Dazap1", "Mus81", "Rps15a", "Sec23b", "Srsf1", "Ppp6r2", "Bsg", "Gm40812", "Med9", "Pgam1", "Irak1", "Gm3604", "Hnrnpa2b1", "Slpi", "LOC108167334", "Dmtf1", "Mapk1", "LOC102633308", "Glrx3", "Hmgn1", "Tmem252", "Ogt", "Man2b1", "Ccnl1", "Pycard", "Racgap1", "Mtdh", "Cox19", "Klhdc4", "Mitd1", "Gart", "Eif3g", "Sfpq", "Hmgb2", "Mgea5", "Gmip", "Clptm1l", "Slc8b1", "Mrps17", "Atxn2l", "Cdca5", "Gm41912", "LOC102632541", "Mterf3", 
## "Nfe2l2", "Rap1a", "Gnas", "0610009B22Rik", "Dync1h1", "Ddb1", "Gm41027", "Acer3", "Eif3a", "Elovl5", "Gm33100", "Grb2", "Slc8a1", "5730408K05Rik", "Sumo3", "Fnip1", "Psme3", "Gm40372", "Zkscan3", "Eif1", "Gm6158", "Ythdf2", "Ms4a2", "Hsp90b1", "Dck", "Ganab", "Arhgap15", "Sdhd", "Zdhhc21", "Aggf1", "Gm40457", "Sec61a1", "Cbx4", "Etfb", "Sltm", "Dek", "Gm38483", "Eif4a1", "Haus4", "Mrpl30", "Gm32137", "Tchh", "Hist1h2br", "Jmjd8", "Ap3s1", "Mapk1ip1l", "Snrnp70", "Pdcd10", "Gm39850", "Cfap36", "Mettl5", 
## "Mrto4", "Ube2c", "Pcgf2", "Mxra7", "Tmco1", "Rara", "Msl1", "Zfp36", "Psen1", "Rbm3", "Nt5c3b", "Mlf2", "S100a13", "Crk", "Zfp719", "Vmp1", "Arl6ip6", "Cd1d2", "Bloc1s2", "Crlf3", "Iqgap2", "Gm30073", "Eid1", "Cdc5l", "Hook3", "Prim1", "Mcm4", "Cnot2", "Minos1", "Tomm5", "Iars", "Pnisr", "Fbxl12", "Rpn1", "Atp6ap2", "Gm39871", "LOC105246245", "Rnf7", "Rp2", "Sdf2", "Ssu72", "Gm32549", "B930092H01Rik", "Glod4", "Osbpl8", "Cfap20", "Nup160", "Tcp10c", "Anxa11", "Bckdhb", "Gadd45b", "Sptssa", "Zfp607", 
## "P4hb", "Ost4", "Gaa", "Poll", "Vmn2r45", "Lman2", "Pes1", "Hyi", "Capn7", "Serinc1", "Eif1ad", "Ilf3", "Atxn7l3b", "Mrps24", "Ndel1", "Oser1", "Plac8", "Rab1a", "Zcchc9", "Ddx21", "Gm38808", "Ywhae", "Gm13429", "Hist1h2ap", "Tmem19", "LOC102633666", "Yeats4", "Layn", "Gm34292", "Taz", "Pnn", "Gm12132", "Gm38416", "Nus1", "Pkn3", "Ier3ip1", "Ctnnb1", "Scamp2", "Ccdc59", "Bop1", "Gm27003", "Timm22", "Ndufs3", "Dcx", "Gm10184", "Sec22b", "Aldoa", "Slc35a4", "Rpl36a", "Akap5", "Catsper4", "Dennd4a", 
## "Eif4ebp2", "Gm37053", "Gm8430", "Mar-06", "Mpp2", "Rpl38", "Tcp1", "Eif2s3x", "Gm38621", "Camk1d", "D17H6S56E-5", "Bnip3l-ps", "Sat1", "Ankrd17", "Dip2b", "Ngdn", "Srp9", "Pak1ip1", "BC031181", "Gpbp1l1", "H2afx", "LOC108168239", "Rps19", "Coro2a", "Son", "Krr1", "Med8", "Rrp9", "Tm2d2", "Naa20", "Paox", "Ehmt1", "Papola", "Slc35c1", "Gm33285", "Rab27a", "Rapgef6", "Patl1", "Acly", "Psma3", "Fen1", "Rin2", "Nfe2", "Lhx1", "Psmc3", "Scarb2", "BC037034", "Snrnp27", "Rnasel", "Atp5f1", "Plekhm2", "Tfpi", 
## "Bag6", "Cd69", "Nif3l1", "Pigm", "

## Warning in AddModuleScore(MKP3HTO, features = list(genes_pd_vs_cntrl_pos), :
## Could not find enough features in the object from the following feature lists:
## Attempting to match case...

## Warning in grep(pattern = paste0("^", s, "$"), x = match, ignore.case = TRUE, :
## argument 'pattern' has length > 1 and only the first element will be used

## Warning in grep(pattern = paste0("^", s, "$"), x = match, ignore.case = TRUE, :
## argument 'pattern' has length > 1 and only the first element will be used

## Warning in grep(pattern = paste0("^", s, "$"), x = match, ignore.case = TRUE, :
## argument 'pattern' has length > 1 and only the first element will be used

## Warning in grep(pattern = paste0("^", s, "$"), x = match, ignore.case = TRUE, :
## argument 'pattern' has length > 1 and only the first element will be used

## Warning in grep(pattern = paste0("^", s, "$"), x = match, ignore.case = TRUE, :
## argument 'pattern' has length > 1 and only the first element will be used

## Warning in grep(pattern = paste0("^", s, "$"), x = match, ignore.case = TRUE, :
## argument 'pattern' has length > 1 and only the first element will be used

## Warning in grep(pattern = paste0("^", s, "$"), x = match, ignore.case = TRUE, :
## argument 'pattern' has length > 1 and only the first element will be used

## Warning in grep(pattern = paste0("^", s, "$"), x = match, ignore.case = TRUE, :
## argument 'pattern' has length > 1 and only the first element will be used

## Warning in grep(pattern = paste0("^", s, "$"), x = match, ignore.case = TRUE, :
## argument 'pattern' has length > 1 and only the first element will be used

## Warning in grep(pattern = paste0("^", s, "$"), x = match, ignore.case = TRUE, :
## argument 'pattern' has length > 1 and only the first element will be used

print(VlnPlot(MKP3HTO, features = "genes_pd_vs_cntrl_pos1", group.by ="CSclassification"))

print(VlnPlot(MKP3HTO, features = "genes_pd_vs_cntrl_pos1", split.by ="CSclassification"))

print(VlnPlot(MKP3HTO, features = "genes_pd_vs_cntrl_pos1", group.by ="CSclassification", split.by = "Idents"))

print(VlnPlot(MKP3HTO, features = "genes_pd_vs_cntrl_pos1", group.by = "Idents"))

projecting the common genes between bulk sequencing DE genes upregulated for PD vs CNTRL and downregulated for PDC+PD VS CNTRL onto the single-cell clusters using the AddModuleScore function and visualised the results using violin plots

genes_pd_vs_cntrl  <- subset(bulk_de_genes, log2FC_PD_vs_Ctrl > 0)
genes_pd_vs_cntrl = rownames(genes_pd_vs_cntrl)
genes_pdpdc_vs_pd = subset(bulk_de_genes, log2FC_PD_DT_vs_PD < 0)
genes_pdpdc_vs_pd = rownames(genes_pdpdc_vs_pd)

final_genes = intersect(genes_pd_vs_cntrl, genes_pdpdc_vs_pd)

MKP3HTO <- AddModuleScore(MKP3HTO,
                    features = list(final_genes), name = "final_genes")

## Warning: The following features are not present in the object: 1810021M19Rik,
## D2Wsu81e, Gm40508, Shfm1, LOC108167358, C330027C09Rik, Gm15264, Sgol1,
## Mir703, Gm40227, Vimp, Gm5812, Gm4617, Tdpx-ps1, 2700094K13Rik, Gm36368,
## Helt, LOC105247125, Gm42312, 5730422E09Rik, Gm36876, Yy2, Gm31619, Ddx26b,
## Gm41847, Gm32673, Fbxl16, LOC108168971, Fdx1l, Rplp0, Gm39923, Gm13394, Tusc5,
## Gm15452, Gm40017, Gm3740, LOC108169013, Gm36673, Rps15a-ps4, Gm14176, Gm6361,
## LOC108167318, 1110018N20Rik, 1810022K09Rik, Rpl37a, Rps2, Gm6532, Gm10705, Rpl8,
## Rps6, Gm39351, LOC108168421, Ube2u, Hn1, Gm32534, Rpl7a, Gm10444, Rps11-ps2,
## Gm32399, Irx5, 1110001J03Rik, Gm39877, Gm30718, Gm11914, Gdap10, Rps12-ps10,
## Gm30181, Gm5620, Rpl13a, Gm41320, Gm3386, Olfr54, 5031434C07Rik, Rps15a,
## Gm40812, LOC108167334, LOC102633308, Tmem252, Mgea5, Gm41912, LOC102632541,
## Gm41027, Gm33100, 5730408K05Rik, Gm40372, Gm6158, Gm38483, Gm32137, Gm39850,
## Cd1d2, Gm30073, Minos1, Gm39871, LOC105246245, Gm32549, B930092H01Rik,
## Tcp10c, Zfp607, Vmn2r45, Hyi, Gm38808, LOC102633666, Gm34292, Gm38416, Dcx,
## Rpl36a, Gm37053, Gm8430, Mar-06, Rpl38, Gm38621, D17H6S56E-5, Bnip3l-ps,
## LOC108168239, Rps19, Gm33285, BC037034, Ppy, Gm6644, Gm29879, Gm8784, Gm34496,
## Gm32269, Gm42042, Gm12911, Gm31727, Gm38663, Gm18798, BE949265, Optc, Tceb1,
## 4931431B13Rik, Gm12338, Gm8995, LOC108167687, Gm7669, LOC105244208, Fam46a,
## Cyp2a22, Gm38761, Rpl18, Gm8866, BC002163, Rps4x-ps, Fam64a, D130037M23Rik,
## Rpl4, Gm5087, Xrcc6bp1, Gm35321, 6330407A03Rik, Gm41192, C330006A16Rik,
## Gm5347, Gm35449, LOC108168951, Gm31842, Peo1, Fam131c, Rpl17, Gm33097,
## Gm29826, Bzrap1, Wbscr22, 2810474O19Rik, Gm33046, Tmem194, Scrt1, Enthd2,
## Emc8-1190005i06rik, Gm10231, Rps3, Gm17057, Vmn1r-ps76, Gm31721, Chrna2,
## Whsc1, Gm29824, 4930543N07Rik, Gm31012, Gm3776, Gm10033, Gm26849, Gm41527,
## Rpl35, Rps18, Gm30145, Gm36754, Gm21158, Gm8463, Gm33651, LOC108167923, Rps15,
## Tbx20, Gm41798, Gm33055, Gm29938, Gm38960, LOC108167404, LOC108168071, Gm20850,
## LOC102634904, Gm32714, Gm42110, Gm31724, Vmn1r-ps78, Gm39697, 0610037L13Rik,
## Rpl23, Gm38957, Gm31208, Fam208b, Gm26688, Gm40403, Fam188a, 1500011K16Rik,
## Smek1, Gm31217, Gm32957, Rps16-ps2, Gm2163, LOC105244649, Sssca1, 2700060E02Rik,
## Gm31758, Obox3-ps5, Hpd, LOC108168233, Gm9892, Gm17250, Gm4332, LOC108167728,
## Gm30346, Gm15710, Rpl13, Gm4754, Gm14680, Gtsf1l, Hmha1, Gm38450, Gm3724,
## Gm39827, Cyp2b26-ps, Gm5486, Lexm, Wfdc6a, Gm5130, LOC108168025, Gm41549,
## Gm18097, Gm5879, Gm4792, Gm9794, Pax4, Gm32394, Cdrt4, Gm4149, Sprr2a2,
## Sepn1, Gm30789, Gm30532, Gm42201, Rpl5, Nodal, Vmn2r44, Gm41519, Rps15a-ps6,
## Rtbdn, Gm4739, Gm39527, LOC102638891, Hbb-bh2, LOC101056014, Gm42213, Gm29933,
## Gm7536, Rps29, Trmt112-ps2, Gm42034, Trim50, LOC102632031, Gm14303, Gm10774,
## Gm39162, LOC108169055, Gm41569, Kirrel2, Serpinb9c, Gm31228, LOC108167924,
## Gm35644, Tmem55b, Gm19600, Gm10222, Gm7180, 2310036O22Rik, Gm14204, Gm9255,
## Rpl34, LOC108167523, Rpl10, Mum1, Gm41220, Ube2n-ps1, Gm3047, Rps26,
## LOC108167669, Fhad1os2, Gm38566, Gm42324, LOC108168511, Gm42344, Gm30330,
## LOC105246089, Gm32727, Gm9745, Rpl39, Rps27a, Olfr287, AV320801, Pnma5,
## Gm42308, 4833422M21Rik, Gm42013, Npm3-ps1, Gm32190, Atp5sl, C030037F17Rik,
## Gm38636, Gm5781, Gm16350, Rps15a-ps7, Gm35343, Pla2g2a, Rps21, Gm16365, Gm32709,
## Gm35396, Gm32756, LOC102640520, Rpl12, Gm33083, Tceb2, Gm33733, LOC108167852,
## Gm39586, Rpl3, Gm6485, 1810043H04Rik, 1700024F13Rik, Gm32670, Gm44504, Gm6988,
## Cyp2c52-ps, Gm38599, Ankrd34a, Gm34521, Gm41302, Elfn2, Rpl14, Gm31048, Barx1,
## Snap25, 9030622O22Rik, Gm33433, 2610524H06Rik, LOC100503338, Gm41204, Gm15500,
## LOC108167658, Selk, Gm30461, Gm41313, Gm42060, Mettl20, Gm34689, Pyhin1,
## Scgb2b8-ps, Gm30658, Rpl19, Gm41166, Gm39665, Igkv3-12, Gm30970, Gm40009,
## Rps10, Fam63a, 5730488B01Rik, Gm34995, Gm35365, Gm16432, Gm40246, Gm13142,
## Rpl37rt, Gm40278, 3110002H16Rik, Gm41995, Fam175b, Gm15723, LOC108167337,
## Gm35884, Obox3-ps6, Gm30005, Gm5278, Gm39487, LOC108168399, Rps20, Gm38732,
## Gm29688, Gm34350, Trpa1, Gbp11, Gm10619, Gm33624, Myom2, Gm31340, LOC108168109,
## Nhp2l1, LOC108168293, C530044C16Rik, Gm41526, Gm41865, Ang6, Serpina1a,
## Gbp2-ps, Gm32099, Gm32422, Rps23-ps1, Gm19303, Gm36188, Rps6ka3, Gm30724,
## Gm40034, Gm7730, 2810025M15Rik, Fam65a, Gm17733, Gm5740, Aes, Gm42398, Gm34865,
## LOC102640586, 1700034E13Rik, Gm2773, Gm41336, Gm13420, Gm8451, Rps9, Gm30364,
## 1810026J23Rik, Gm35208, LOC100534331, LOC108169173, Gm11478, Rps14, Gm40799,
## Gm30286, Gm39866, Gm14199, Gm32525, Gm19764, Morf4l1-ps1, Gm38782, Gm7068,
## Vmn1r-ps71, Gm7666, Gm12371, Gm35040, LOC108168082, Dpt, Gm9703, Ugt2b36,
## D17Wsu92e, Gm30212, Gm40191, Rpl9-ps6, Gm40446, Rps12-ps4, Gm40680, Zufsp,
## LOC108167317, Rps13-ps2, Gm8159, D230017M19Rik, Amd-ps4, Zic3, Gm14586, Rpl11,
## Gm20687, 9330133O14Rik, Gltscr1l, Gm32335, Gm15163, Tcp10a, Snhg7, Gm32370,
## Gm4855, Gm32896, Gm38554, Gm4956, Gm35760, Gm28166, Gm32276, Kdelc1, Gm20836,
## Cyp2j11, Gm32450, Rplp2, Gm34070, Gm9843, Gm34472, LOC108167644, LOC102637814,
## Gm1966, Gm33531, Gm30628, Gm2999, Gm36670, Rpl41, Gm3086, Gm14494, Gm40733,
## Gm30524, Gm14532, Gm11683, Cecr5, Gltscr2, Gm31546, LOC108168860, Gm39101,
## Gm30215, Myocd, Gm35118, Pcsk2os1, Gm15610, Gm4827, Ngfrap1, Gm31102, Gm32545,
## Gm35078, Gm40905, LOC108168267, Rps15a-ps5, Gm30548, LOC108169051, Rps27a-ps2,
## Gm27177, Rps12, Gm30368, Csrp2bp, Cyp2j15-ps, Rpl21, Rps16, Casc5, Gm18075,
## Gm16701, Prap1, Atp4b, Gm40874, Rn18s-rs5, Csta1, Rpl10a-ps1, Fgf4, Pcdhga3,
## LOC105245440, Gm34317, Gm8267, Gm30400, Gm33983, Pisd-ps2, 4930432L08Rik,
## 9430014N10Rik, Gm36603, Vmn1r63, LOC108167528, Gm42088, LOC108169076,
## Gm39160, Krt88, Cfap161, Gm35478, Rplp1, Rpl30, Gm39325, Gm12496, Gm13268,
## Kcnq4, Kdelc2, Btnl6, Gm3219, Gm11682, Gm32017, LOC108167997, Gm39481, Pex5l,
## Suv420h1, LOC107988026, Morf4l1b, Fam109b, Gm16348, 6820431F20Rik, Gm34333,
## LOC105246961, Asun, LOC102641351, Gm33893, Hottip, Gm10123, Gm33068, Gm8926,
## Gm33950, Gm36025, Prl7a1, LOC108169047, A430046D13Rik, Gm10857, Zscan30,
## Gm42235, Gm32110, Rpl7, Igf2bp1, Gm31042, Nat6, D930030F02Rik, LOC105243337,
## Gm2803, Gm41508, Gm9712, Fam60a, Tpsab1, LOC102634709, LOC108168357, Gm13528,
## Gm31962, Rpl18-ps1, Gm41401, 6330416G13Rik, Platr11, 4833447I15Rik, Gm32348,
## Gm42259, Gm5907, Gm9769, LOC108168928, Rpl7l1, Gm36048, Gm40336, Gm2477,
## LOC102631912, Gm38512, Gm39033, Fam134a, Gm36220, 4933427G17Rik, Gm35104,
## Gm41892, LOC108167447, LOC100534390, LOC105244803, LOC108167686, 1700112E06Rik,
## Gm40363, 6430584L05Rik, Gm31078, Gm41209, Mar-08, Synj2bp-cox16, Gm30140,
## Gm40159, Gm26232, Olfr1418, Gm13248, Dgcr14, Gm16505, Gm34388, Gm38515,
## 5930403L14Rik, Gm6498, Gm39851, Fam179b, Gm36596, LOC108168889, 4933413G19Rik,
## Cmtm5, Dirc2, Gm34058, Gm41926, Gm34667, Gm30352, Gm38767, Gm15387, Myh7,
## Ang-ps2, 4831440E17Rik, Wnt2, Gm35765, Gm34444, Gm12119, 4931406H21Rik, Gm40349,
## LOC108167465, 1700112H15Rik, Pllp, 1700026F02Rik, Rps17, Gm10221, Gm35155,
## Prss43, Adh4, Gm32374, Nudc-ps1, Sox10, Gm39974, Gm4184, 1700031M16Rik, Gm39337,
## LOC108167633, Gm2048, Gm5682, Etd, Kif26a, Gm32579, LOC102642832, 1700018B08Rik,
## LOC108167604, Fgb, Mgat5b, Gm32328, Gm7426, Gm9835, LOC108168813, LOC108167469,
## Gm42325, Gm18837, LOC105242790, Hoxc5, LOC108167878, Itln1, Gm32220, Rpl6,
## Opalin, Fam96b, Nup62cl, C330016L05Rik, Gm15742, Gm31244, Gm31473, Gm42250,
## Ugt1a10, Gm38499, Rpl31, LOC105243608, Cacna1h, Gm6650, Gm38993, Gm34669,
## Best2, Gm40550, 1110008F13Rik, Rpl10-ps3, Gm31298, Hils1, C230062I16Rik,
## Gm30936, Mir6236, Gm29955, Tekt3, Gm20684, LOC108168493, Rps11, Gm30086,
## Gm39241, Mar-05, LOC108169093, Speer4a, Gm7265, Stoml3, Gm41836, 3110062M04Rik,
## Tmem114, LOC102632770, Gm7338, LOC102640103, Gm41359, Gm4726, Rps3a1,
## Gm12346, Gm31265, Gm32646, Olfr129, Ssfa2, Serpinb1b, A730089K16Rik, Gm41699,
## Gm9207, 4930447F24Rik, Gm4750, Gm10224, LOC108168488, Gm36386, LOC108167736,
## LOC108167896, Pcp4, Rpl7a-ps12, Wnt1, Defb26, Gm15612, Gm30922, LOC100503946,
## Mtl5, Gm3050, LOC102642435, Gm29943, Spint3, LOC102635252, Rps8, Gm34016,
## Gm36229, Gm18221, Rpl28-ps1, Gm9919, Gm39938, Gm39546, Rps27rt, Gm32069,
## Gm30848, LOC108167773, Gm41598, Pcdha11, 1700003C15Rik, Fam65b, Gm34749,
## Gm40986, 6030466F02Rik, Gm33236, Gm41056, Grxcr2, Amica1, Fam69a, Btnl7-ps,
## Mtss1l, Gm31671, Olfr671, Sep-15, 4930447M23Rik, Rpsa, Raver1, 4833427

MKP3HTO = subset(x = MKP3HTO, subset = CSclassification=="pDC-depletion", invert = TRUE)
Idents(MKP3HTO) = MKP3HTO$Idents
print(VlnPlot(MKP3HTO, features = "final_genes1", group.by ="CSclassification"))

print(VlnPlot(MKP3HTO, features = "final_genes1", split.by ="CSclassification"))

print(VlnPlot(MKP3HTO, features = "final_genes1", group.by ="CSclassification", split.by = "Idents"))

print(VlnPlot(MKP3HTO, features = "final_genes1", group.by = "Idents"))

p = VlnPlot(MKP3HTO, features = "final_genes1", group.by = "Idents")
print(p)

#ggsave("~/Downloads/final_genesnew.svg", plot = p, width = 18, height = 10, limitsize = FALSE)

This block of code finds the differentially expressed genes for each cluster between the three conditions (Bl6 vs PD, PD vs pDC-PD, and Bl6 vs pDC-PD) using the wilcox test for all the genes of the seurat object. The DE genes are selected based on the threshold of avg_log2FC > 0.25 and p_val_adj < 0.05. The DE genes are then plotted onto a volcano plot using the package EnhancedVolcano. Finally, the DE genes are run through enrichR for a gene set enrichment analysis.The genes are filtered to a positive set if their log fc > 0 and a negative set if their log fc < 0 for separating the input to the enrichment analysis.

Idents(MKP3HTO) = MKP3HTO$CSclassification

# Define the conditions and cluster names
conditions <- c("Bl6", "Platelet-depletion", "pDC-Platelet-depletion")
cluster_names <- c("0", "1", "2", "3", "4", "5", "6", "7", "8", "9")

# Perform differential expression and enrichment analysis for each cluster
for (cluster_name in cluster_names) {
  cat("Cluster:", cluster_name, "\n")
  
  # Define a list of comparisons for each cluster
  comparisons <- list(
    "Bl6_vs_Platelet-depletion" = c("Bl6", "Platelet-depletion"),
    "Platelet-depletion_vs_pDC-Platelet-depletion" = c("Platelet-depletion", "pDC-Platelet-depletion"),
    "Bl6_vs_pDC-Platelet-depletion" = c("Bl6", "pDC-Platelet-depletion")
  )
  subset_obj = subset(x = MKP3HTO, subset = seurat_clusters == cluster_name)
  for (comparison_name in names(comparisons)) {
    conditions_to_compare <- comparisons[[comparison_name]]
    
    # Extract gene names
    gene_names <- rownames(subset_obj@assays$RNA@scale.data)[subset_obj$seurat_clusters == cluster_name]
    
    # Perform differential expression analysis
    de_results <- FindMarkers(
      subset_obj,
      ident.1 = conditions_to_compare[1],
      ident.2 = conditions_to_compare[2],
      test.use = "wilcox",
      genes.use = gene_names  # Specify the gene names for the comparison
    )
    
    
    # Get the list of differentially expressed genes
    de_genes <- de_results[abs(de_results$avg_log2FC) > 0.25 & de_results$p_val_adj < 0.05, ]
    # Define the file name
    file_name <- paste("cluster", cluster_name, comparison_name, ".csv", sep = "_")
    de_genes$genes = rownames(de_genes)
    # Write the results to a CSV file
    #write.csv(de_genes, file = file_name, row.names = TRUE)
    
    #cat("Saved results to file:", file_name, "\n")
    
    enrich_genes_neg = rownames(de_genes[de_genes$avg_log2FC < 0, ])
    enrich_genes_pos = rownames(de_genes[de_genes$avg_log2FC > 0, ])
    
    enrichment_results_neg <- enrichR::enrichr(
      genes = enrich_genes_neg,
      databases = "GO_Biological_Process_2023"
    )
    enrichment_results_neg = as.data.frame(enrichment_results_neg)
    enrichment_data_neg <- data.frame(
  Term = enrichment_results_neg$GO_Biological_Process_2023.Term,
  Overlap = enrichment_results_neg$GO_Biological_Process_2023.Overlap,
  P.value = enrichment_results_neg$GO_Biological_Process_2023.P.value,
  Adjusted.P.value = enrichment_results_neg$GO_Biological_Process_2023.Adjusted.P.value,
  Odds.Ratio = enrichment_results_neg$GO_Biological_Process_2023.Odds.Ratio,
  Combined.Score = enrichment_results_neg$GO_Biological_Process_2023.Combined.Score,
  Genes = enrichment_results_neg$GO_Biological_Process_2023.Genes
)
    if (nrow(enrichment_data_neg) == 0) {
          next
    }else{
       enrich_neg = plotEnrich(enrichment_data_neg)
       print(enrich_neg)
      #ggsave(paste("cluster", cluster_name, comparison_name, "neg.png", sep = "_"), plot = enrich_neg)
}
    
    enrichment_results_pos <- enrichR::enrichr(
      genes = enrich_genes_pos,
      databases = "GO_Biological_Process_2023"
    )
    enrichment_results_pos = as.data.frame(enrichment_results_pos)
        enrichment_data_pos <- data.frame(
  Term = enrichment_results_pos$GO_Biological_Process_2023.Term,
  Overlap = enrichment_results_pos$GO_Biological_Process_2023.Overlap,
  P.value = enrichment_results_pos$GO_Biological_Process_2023.P.value,
  Adjusted.P.value = enrichment_results_pos$GO_Biological_Process_2023.Adjusted.P.value,
  Odds.Ratio = enrichment_results_pos$GO_Biological_Process_2023.Odds.Ratio,
  Combined.Score = enrichment_results_pos$GO_Biological_Process_2023.Combined.Score,
  Genes = enrichment_results_pos$GO_Biological_Process_2023.Genes
)
# Plot the enrichment results
        if (nrow(enrichment_data_pos) == 0) {
          next
        }else{
          enrich_pos = plotEnrich(enrichment_data_pos)
          print(enrich_pos)
          #ggsave(paste("cluster", cluster_name, comparison_name, "pos.png", sep = "_"), plot = enrich_pos)
  }
}
}

## Cluster: 0 
## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.

## No genes have been given

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.
## Cluster: 1

## No genes have been given

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.

## No genes have been given

## Cluster: 2 
## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.

## Warning in plotEnrich(enrichment_data_pos): There are duplicated trimmed names
## in the plot, consider increasing the 'numChar' setting.

## No genes have been given

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.

## Warning in plotEnrich(enrichment_data_pos): There are duplicated trimmed names
## in the plot, consider increasing the 'numChar' setting.

## Cluster: 3 
## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.

## Warning in plotEnrich(enrichment_data_neg): There are duplicated trimmed names
## in the plot, consider increasing the 'numChar' setting.

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.

## Warning in plotEnrich(enrichment_data_neg): There are duplicated trimmed names
## in the plot, consider increasing the 'numChar' setting.

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.

## Cluster: 4 
## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.
## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.
## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.

## Cluster: 5 
## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.

## No genes have been given
## No genes have been given

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.
## Cluster: 6

## No genes have been given
## No genes have been given

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.

## No genes have been given

## Cluster: 7

## No genes have been given

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
## Parsing results... Done.

## Cluster: 8

## No genes have been given
## No genes have been given
## No genes have been given

## Cluster: 9

## No genes have been given
## No genes have been given
## No genes have been given

scatterplots to plot the common DE genes between single-cell clusters (MK-MEP cycling and MK-MEP) and bulk-seq for the comparison PD vs CNTRL. the axis’ represent the expression levels of the genes in the experiments.

cluster_3_Platelet_depletion_vs_Bl6_common_genes_pos = read.csv("~/Downloads/cluster_3_Platelet-depletion_vs_Bl6_common_genes_pos.csv")
cluster_3_Platelet_depletion_vs_Bl6_common_genes_pos = cluster_3_Platelet_depletion_vs_Bl6_common_genes_pos$x


subset_obj = subset(x = MKP3HTO, subset = seurat_clusters == 3)
Idents(subset_obj) = subset_obj$CSclassification
de_results <- FindMarkers(
      subset_obj,
      ident.1 = "Platelet-depletion",
      ident.2 = "Bl6",
      test.use = "wilcox",
      genes.use = cluster_3_Platelet_depletion_vs_Bl6_common_genes_pos  # Specify the gene names for the comparison
    )

de_results = de_results[cluster_3_Platelet_depletion_vs_Bl6_common_genes_pos, ]    
genes_pd_vs_cntrl_pos  <- subset(bulk_de_genes, log2FC_PD_vs_Ctrl > 0.25 & (FDR_PD_vs_control) < 0.05)
genes_pd_vs_cntrl_pos = genes_pd_vs_cntrl_pos[, c("log2FC_PD_vs_Ctrl", "FDR_PD_vs_control")]
genes_pd_vs_cntrl_pos = genes_pd_vs_cntrl_pos[cluster_3_Platelet_depletion_vs_Bl6_common_genes_pos, ]
VlnPlot(subset_obj, features = cluster_3_Platelet_depletion_vs_Bl6_common_genes_pos)

genes_pd_vs_cntrl_pos$Gene = rownames(genes_pd_vs_cntrl_pos)
de_results$Gene = rownames(de_results)

merged_data <- merge(de_results, genes_pd_vs_cntrl_pos, by = "Gene", all = TRUE)
rownames(merged_data) = merged_data$Gene
# Create the scatterplot
p = ggplot(data = merged_data, aes(x = avg_log2FC, y = log2FC_PD_vs_Ctrl, label = Gene,color = avg_log2FC)) +
  geom_point() +
  geom_text(hjust = 0.001, vjust = 0.001, nudge_x = 0.001, nudge_y = 0.001, size = 3) + # Adjust label position and size
  scale_color_gradient(low = "blue", high = "red") +
  labs(x = "Single Cell Cluster MK-MEP_cycling", y = "Bulk Seq PD vs Bl6") +
  ggtitle("Scatterplot of Fold Changes for common PD vs Bl6 DE genes")
print(p)

#ggsave("scatter_plot1.svg", plot = p, width = 15, height = 8, limitsize = FALSE)

## response to type I interferon genes as bulk_seq genes
interferon = read.csv("~/Downloads/GOBP_RESPONSE_TO_TYPE_I_INTERFERON.csv", header = TRUE)


## converting the genes to mouse gene form - "Tyk2" "Hla-a"
interferon$PROBE = tolower(interferon$PROBE)
firstup = function(x) {
  substr(x, 1, 1) = toupper(substr(x, 1, 1))
  
  x
}
interferon$PROBE = firstup(interferon$PROBE)
rownames(interferon) = interferon$PROBE
interferon = head(rownames(interferon), 40)

genes_pd_vs_cntrl_pos  <- bulk_de_genes
genes_pd_vs_cntrl_pos = genes_pd_vs_cntrl_pos[, c("log2FC_PD_vs_Ctrl", "FDR_PD_vs_control")]
genes_pd_vs_cntrl_pos = genes_pd_vs_cntrl_pos[interferon, ]
genes_pd_vs_cntrl_pos = na.omit(genes_pd_vs_cntrl_pos)
inf_final = rownames(genes_pd_vs_cntrl_pos)
genes_pd_vs_cntrl_pos$Gene = rownames(genes_pd_vs_cntrl_pos)


# Perform differential expression analysis using the binary variable
subset_obj = subset(x = MKP3HTO, subset = seurat_clusters == 3)
Idents(subset_obj) = subset_obj$CSclassification
de_results_inf <- FindMarkers(subset_obj, ident.1 = "Platelet-depletion", ident.2 = "Bl6", test.use = "DESeq2")

## converting counts to integer mode

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

de_results_inf = de_results_inf[inf_final, ]
de_results_inf = na.omit(de_results_inf)
de_results_inf$Gene = rownames(de_results_inf)


merged_data <- merge(de_results_inf, genes_pd_vs_cntrl_pos, by = "Gene", all = TRUE)
rownames(merged_data) = merged_data$Gene
# Create the scatterplot
p = ggplot(data = merged_data, aes(x = avg_log2FC, y = log2FC_PD_vs_Ctrl, label = Gene,color = avg_log2FC)) +
  geom_point() +
  geom_text(hjust = 0.001, vjust = 0.001, nudge_x = 0.001, nudge_y = 0.001, size = 3) + # Adjust label position and size
  scale_color_gradient(low = "blue", high = "red") +
  labs(x = "Single Cell Cluster MK-MEP_cycling", y = "Bulk Seq PD vs Bl6") +
  ggtitle("Scatterplot of Fold Changes for interferon response genes PD vs Bl6")
print(p)

## Warning: Removed 5 rows containing missing values (`geom_point()`).

## Warning: Removed 5 rows containing missing values (`geom_text()`).

#ggsave("scatter_plot3.svg", plot = p, width = 20, height = 8, limitsize = FALSE)

dotplot of genes per category to determine the clusters more accurately

# Load the required Seurat library


# Load your Seurat object (replace 'seurat_obj' with your Seurat object)
seurat_obj <- MKP3HTO

# Define the gene sets as a list of lists
    MK_TF = c("Cited2", "Fli1", "Pbx1", "Mef2c", "Meis1")
    MK_GRANULES = c("Pf4", "Vwf", "Clu", "Ppbp", "Serpine2", "Thbs1", "F5", "Ctla2a", "Angpt1")
    MK_SIGNALING = c("Rap1b", "Cavin2", "Tpm4", "Myl9", "Gng11", "Alox12", "Tmsb4x", "Plek", "Prkg1", "Prkca", "Trpc6")
    MK_SURFACE = c("Itga2b", "Mpl", "Cd9", "Gp9", "Gp1bb", "Gp5", "Gp6", "Clec1b", "Slamf1")
    ERY_TF = c("Zfpm1", "Hmbg3", "Gata1", "Gfi1b", "Klf1")
    ERY = c("Eng", "Car2", "Hba-a2", "Gypa")
    GMP_TF = c("Gata2", "Cebpa", "Irf8")
    GMP = c("Elane", "Mpo")
    NEUTRO_TF = c("Spi1", "Cebpe")
    NEUTRO = c("Camp", "Ltf", "Retnlg")
    NK_TF = c("Ets1", "Foxo1", "Runx3")
    NK = c("Skap1", "Ccl5", "Id2")
    BASO_TF = c("Lmo4", "Runx1", "Stat5b")
    BASO = c("Hdc", "Cd200r3", "Prss34")
    Cell.cycle = c("Nusap1", "Top2a", "Mki67")
    MIXED_RB = c("Ncl", "Npm1", "Cmss1")
    MIXED_TI = c("Eif4b", "Eif4a1", "Eif2s2", "Pdcd4")
    MIXED_NM = c("Nme1", "Dctpp1")
    MIXED_PC = c("Hsp90aa1", "Hsp90ab1")

features <- list("MK TF" = MK_TF, "MK_GRANULES" = MK_GRANULES, "MK_SIGNALING" = MK_SIGNALING, "MK_SURFACE" = MK_SURFACE, "ERY_TF" = ERY_TF, "GMP_TF" = GMP_TF,"GMP"=GMP,"NEUTRO_TF" = NEUTRO_TF,"NEUTRO" = NEUTRO, "NK_TF"=NK_TF,"NK" = NK,"BASO_TF" = BASO_TF,"BASO" = BASO, "Cell.cycle" = Cell.cycle,"MIXED_RB"=MIXED_RB,"MIXED_TI" = MIXED_TI, "MIXED_NM" = MIXED_NM, "MIXED_PC" = MIXED_PC)

d = DotPlot(object = MKP3HTO, features=features, cluster.idents=T) + theme(axis.text.x = element_text(angle = 90))

## Warning in FetchData.Seurat(object = object, vars = features, cells = cells):
## The following requested variables were not found: Hmbg3

## Warning: Scaling data with a low number of groups may produce misleading results

## Warning: The `facets` argument of `facet_grid()` is deprecated as of ggplot2 2.2.0.
## ℹ Please use the `rows` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

print(d)